Semi-analytical design of antireflection gratings for photonic crystals 
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This article concerns the design of antireflection structures which, placed on a photonic crystal 
surface, significantly diminish the fraction of energy lost to reflected waves. After a review of the 
classes of these structures proposed to date, a new method is presented in detail for the design of 
antireflection gratings operating in a wide range of angles of incidence. The proposed algorithm is 
illustrated by means of several examples, showing the advantages and limitations. 



I. INTRODUCTION 

Photonic-crystal (PC) devices may become vital ingre- 
dients of integrated optical circuits. For instance, they 
open new possibilities and functionalities like negative 
refraction 1-3 , supercollimation 4,5 , compact isolators 6 ' 7 , 
et c. On the other hand, PCs generally present significant 
reflectivity at their surface while, in many applications, 
one strives after a perfect energy transfer between the 
incident plane wave and the propagative PC cigenmode 
(or, more rarely, several such modes). An important issue 
is then to eliminate propagative reflected waves. In this 
paper we propose a semi-analytical algorithm for the de- 
sign of gratings acting as wide-angle antireflection (AR) 
structures for two-dimensional PCs. 

The article is organized as follows. In the next section 
we review existing AR structures for homogeneous media 
and PCs. Subsequently, in section III, the three steps of 
the new design procedure are presented in detail. This 
algorithm is then tested on several particular cases: a PC 
flat lens (section IV), a supercollimating PC (section V), 
and a non- reciprocal PC mirror (section VI). Finally, in 
section VII we discuss briefly the influence of AR gratings 
in the evanescent- wave regime. 



II. BACKGROUND 

In this section we shall review existing methods of elim- 
inating reflections from interfaces which separate differ- 
ent media. We start with the simpler case of interfaces 
separating homogeneous media, since the basic tools used 
to deal with this class of systems turn out to be valu- 
able also when more complex media, such as PCs, are 
involved. The existing AR solutions can be divided into 
three broad classes. The refractive index n of AR coatings 
depends only on the coordinate z perpendicular to the in- 
terface; depending on whether n(z) is piecewise-constant 
or not, we speak of homogeneous-layer or inhomogeneous- 
layer AR coatings 8 . In turn, the refractive index of AR 



gratings is also a function of the coordinates parallel to 
the interface. Figure 1 shows example AR structures be- 
longing to these three classes. 

The simplest homogeneous-layer AR coating consists 
of a single thin film with refractive index n and thick- 
ness d chosen so as to ensure destructive interference of 
waves reflected from its top and bottom surfaces, excited 
by a plane wave with a certain frequency lo and angle of 
incidence 9. The values of n and d can be obtained ana- 
lytically. In order to improve the angular and frequency 
tolerance of the coating, one can increase the number of 
layers to make the transition between the refractive index 
of the super- and substrate more gradual 8 , thus reducing 
the amplitude of waves reflected on the individual dis- 
continuities of n(z). In the limit of an infinite number 
of layers, one arrives at an inhomogeneous-layer coating 
with a continuous monotonic profile n(z). Several spe- 
cific profiles have been proposed in the literature (see ref. 
9 for a review). Unfortunately, these "ideal" AR coatings, 
even the single-layer one, cannot usually be realized be- 
cause of the lack of suitable materials with the required 
values of n. This is the case, in particular, for air-solid 
interfaces, where thin films with refractive index close to 
unity are needed. 

Short of using special ultra-low-index media, such as 
so-called Reststrahlen materials 8 , there are two basic 
ways of circumventing this problem. One consists in 
abandoning completely the structures based on a "con- 
tinuous transition" between the super- and substrate re- 
fractive indices in favor of interference-based coatings 8 . 
In these systems, the total reflected wave vanishes thanks 
to the destructive interference of partial waves generated 
at interfaces between layers with contrasting values of n. 
For instance, so-called w-coatings consist of two layers 
with prescribed refractive indices n\ and n 2 (correspond- 
ing to realistic materials) and thicknesses d\ and d 2 ad- 
justed so as to eliminate reflection at the desired values 
of lo and 0. Note that n\ and n 2 must satisfy certain con- 
ditions in order that appropriate d\ and d 2 can be found 
(ref. 10, p. 96). A disadvantage of v-type AR coatings is 
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FIG. 1: Example AR structures belonging to the three principal classes described in the text. Darker areas denote regions 
with higher refractive index, (a) A homogeneous-layer AR coating, (b) An inhomogeneous-layer AR coating, (c) A (binary 
lamellar) AR grating. 



that their angular and frequency tolerance are usually in- 
ferior even to those of the corresponding ideal single-layer 
coatings (ref. 10, p. 97; rcf. 11, p. 188). 

The other solution consists in using subwavelength 
gratings to simulate AR coatings with arbitrary n(z) 
profiles. In many cases, the effective-medium theory of 
gratings can be employed to calculate the grating profile 
mimicking the desired n(z) dependence 9 . Several types 
of gratings, such as the lamellar, trapezoidal, sinusoidal, 
triangular and pyramidal ones, have been studied in the 
literature and shown to have good AR properties 9 ' 12-14 . 
A review of the experimental methods used to fabricate 
such AR structures can be found in ref. 15. 

Let us now turn to the case of PCs. One of their dis- 
tinguishing features is the dependence of their reflection 
coefficient on the position of their truncation plane. One 
could hope then that a significant reduction of a PC's re- 
flectance could be achieved without adding any AR struc- 
ture, but simply by choosing an appropriate cut. For 
some crystals, this has indeed proved to be possible 16 ' 17 . 
In particular, Botten et al. 17 have shown that very low 
reflectance is a rather general feature of rod-type PCs 
truncated midway between successive layers of rods. For 
many crystals, however, no truncation plane provides a 
sufficiently small value of reflectance (see fig. 2). This 
method of reducing reflection is therefore not general 
enough, and one often has to resort to introducing some 
AR structure. Several types of them have been proposed 
in the literature. More often than not, they have much 
in common with one of the solutions developed with ho- 
mogeneous materials in mind, reviewed in the previous 
subsection. 

The simplest approach has been proposed by Li et 
al. 18 . It consists in placing in front of the crystal a v- 
type AR coating composed of two layers made of con- 
stituent materials of the crystal, usually air and a dielec- 
tric. Their thicknesses can be determined analytically 
or graphically as soon as the reflection coefficient of the 
semi-infinite uncoated PC at the selected operation fre- 
quency is known. This frequency is assumed to be low 
enough that only the zeroth diffraction order be prop- 
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FIG. 2: The dependence of the reflectance \r \ 2 of the PC 
shown in the inset, placed in air, on the location of the trun- 
cation plane z = zq. The PC consists of a hexagonal lattice of 
air holes with radius 0.365a, where a is the lattice constant, 
etched in a dielectric matrix with permittivity e = 10.6. The 
impinging wave is normally incident, s-polarized and has fre- 
quency lo = 0.311 x 2-Kc/a. It can be seen that |r | 2 does not 
fall under 0.13 for any truncation plane. 



agative in the layer closest to the PC surface, so that all 
higher orders reflected by the crystal decay before reach- 
ing the boundary between the two layers. If it is not the 
case, the structure may still work, but the optimum layer 
thicknesses will not in general be given by the analytical 
formulas derived from the thin-film theory. The basic 
disadvantage of this otherwise very attractive approach 
is the relatively low angular and frequency tolerance of 
w-type AR coatings. 

Another group has studied a similar approach, in which 
one of the homogeneous layers is replaced by a grating 
of rods or holes 19 ' 20 . The performance of the resulting 
AR gratings was shown to be similar to that of u-type 
coatings 20 . Related structures, albeit with only a single 
degree of freedom (the radius of the outermost rods of a 
rod-type PC), were also analyzed by Jin and He 21 . In 
turn, Zhang and Li 22 proposed a more complex wide- 
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angle AR grating for the interface between air and a PC 
flat lens, whose fabrication would be seriously hindered, 
though, by the presence of extremely thin air slits and 
dielectric veins. 

In contrast to the above works, which concentrated 
on the low-frequency regime, the authors of refs. 23-25 
endeavored to eliminate the reflection from an interface 
between a PC and a semiconductor at a frequency where 
multiple propagative diffraction orders existed in the lat- 
ter medium. They proposed an AR grating consisting of 
elongated drop-shaped air holes, this nonstandard form 
being motivated by the desire of ensuring a gradual tran- 
sition between the two media. In fact, the shape of the 
resulting structure resembles closely the classical trape- 
zoidal AR gratings studied, e.g., by Raguin and Morris 14 . 
The improvement brought about by these gratings has 
subsequently been demonstrated experimentally 26 . Un- 
fortunately, Baba et al. did not provide any analytical 
guidelines regarding the choice of the geometrical pa- 
rameters of the gratings, resorting instead to a time- 
consuming scan of the parameter space in order to find 
the optimum structure 23,24 . 

A very important contribution was made by Lawrence 
et al. 27 ' 28 , who introduced the concept of generalized 
matrix-valued effective immittance S of gratings and 
showed that their effective- medium description can be 
made arbitrarily accurate by allowing S to have suf- 
ficiently large dimensions. Such matrix-valued immit- 
tance cannot serve as a drop-in replacement of the corre- 
sponding scalar quantity in standard formulas derived for 
homogeneous media 28 . Therefore, for instance, analyti- 
cal determination of the optimum parameters of an AR 
structure for a given PC composed of layers described by 
a matrix S is not possible. However, Lawrence et al. 2S 
derived equations similar to the classical Fresnel formu- 
las, but involving generalized immittances, and showed 
that accurate values of the scattering coefficients of typi- 
cal PCs could be obtained already using immittance ma- 
trices truncated to 5 x 5 elements or less. Since numer- 
ical calculations involving such small matrices are very 
fast, it becomes feasible to design optimum AR grat- 
ings by performing a full scan of the available param- 
eter space. Example AR structures presented in ref. 28 
include, for instance, a relatively wide-band AR coating 
of a superprism-type PC at a frequency where multiple 
propagative diffraction orders exist in the adjacent dielec- 
tric. A slight limitation of the approach of Lawrence et 
al. is that the generalized immittance has only been de- 
fined for 2D gratings symmetrical with respect to a two- 
fold rotation axis parallel to the direction of invariance. 
Therefore, it does not cover, for instance, triangular or 
trapezoidal gratings. 

All the AR structures discussed so far are relatively 
compact, with thickness rarely exceeding one or two lat- 
tice constants of the underlying PC. Some authors have 
advocated sacrificing compactness in favor of potentially 
larger frequency and angular tolerance offered by thick 
stacks of gratings, whose geometry changes gradually 



so as to ensure a smooth ("adiabatic") transition of the 
electromagnetic field of the incident plane wave towards 
the Bloch mode of the semi-infinite PC. Several design 
principles for such gratings have been proposed 29,30 . In 
this work, we shall focus on compact AR structures, and 
therefore we omit a detailed discussion of adiabatic AR 
gratings. 



III. DESIGN PROCEDURE 

In this section we shall present an alternative method 
of designing AR gratings for PCs. Compared to the ap- 
proaches reviewed in section II, it has the following distin- 
guishing features. First, instead of performing a poten- 
tially time-consuming global scan of possible geometries, 
we use an effective-medium model of PCs to calculate 
analytically the geometrical parameters of a "tentative" 
AR grating; then, if necessary, we refine these parame- 
ters with a numerical local-minimization algorithm. The 
final result is a trapezoidal AR grating. Second, the AR 
structures obtained in the proposed way have good an- 
gular tolerance, which in some important applications of 
PCs is more significant than the frequency tolerance. For 
instance, the quality of the image produced by a PC flat 
lens with effective refractive index n « — 1 depends in the 
first place on the angular range of incident plane waves 
which arc transmitted through the lens with little or no 
energy loss. On the other hand, frequency tolerance is 
not vital since n can be close to —1 only in a narrow 
frequency band. In our developments we draw on the 
results of Raguin and Morris 14 , who demonstrated that 
triangular and trapezoidal gratings allow a significant re- 
duction of reflectance at interfaces between homogeneous 
media in a wide range of incidence angles. 

The proposed design algorithm consists of three ba- 
sic steps. In step 1 we use the value of the reflection 
coefficient of the PC at a particular frequency and inci- 
dence angle to calculate analytically the parameters of 
a homogeneous AR coating appropriate for the PC. In 
step 2 we convert the coating into a lamellar AR grating 
composed of the same materials as the PC itself. If nec- 
essary, in step 3 we adjust the shape of the grating using 
a numerical optimisation procedure, obtaining finally a 
trapezoidal AR grating. We shall now proceed to detail- 
ing these three constituent steps of the algorithm. They 
are summarised in the flowchart in fig. 3. 

Step 1 We begin by using the classical theory of AR 
coatings to calculate the refractive index and thickness 
of a single-layer coating that should minimize the reflec- 
tion from the PC surface at a fixed angle of incidence 9 
and frequency uj. The value of 9 should lie approximately 
midway the desired angular operation range of the final 
AR structure. Consider the system shown in fig. 4, in 
which a propagative plane wave with frequency lo = cfco 
and wave vector k\ = (k x ,k z i) — (nikosm6,nikocos6) 
impinges from the dielectric I with refractive index n\ on 
the surface of the homogeneous film 2 with thickness cfe 
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FIG. 3: Successive steps of the proposed AR grating design 
algorithm. 



and refractive index n 2 coating the PC 3. Let us assume 
that medium 2 supports exactly one propagativc diffrac- 
tion order and is sufficiently thick for all the evanescent 
orders reflected from the PC surface to vanish at the in- 
terface between media 1 and 2. We can calculate the total 
amplitude of the reflected plane wave, r, by summing up 
all the multiple reflections occurring in the system: 



r = r 12 +ti2 < Sr23$[^(r2i$r 2 3$)" 
_ 7"12 + (£12*21 - r 12 r 21 )r 23 § 2 
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where and 



denote the amplitudes of the waves 
reflected from the interface between media i and j and 



transmitted through it, respectively, and $ 



2 with 



k z2 = {n^kQ — kD 1 / 2 . The coefficients r i:) and tij are given 
by the Fresnel's formulas, 
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with Ei denoting the transverse immitance of medium i 
normalised to the immitance of free space. For s po- 
larisation Si is defined as the transverse impedance of 
medium i: S, = = fc ^i/fc Zii . For p polarisation, it 
is defined as the transverse admittance of the medium: 
Ei = Yi = koCi/k z ,i. From eqs. (2) it immediately follows 
that r 2 i = — ri2 and £12*21 - ri 2 r 2 i = 1, hence 



ri2 + r 23 $ 2 
1 + ri 2 r 23 $ 2 ' 



(3) 



The parameters of the antireflection coating, n 2 and d,2, 
can now be obtained by requiring the numerator of the 
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FIG. 4: System considered in step 1 of the AR grating design 
procedure. 



fraction in the above equation to vanish. If the coating 
is lossless, so that |$| = 1, the numerator vanishes if and 
only if (i) the moduli of r i2 and r23 are equal and (ii) the 
thickness d 2 is such that 



argr 2 3 + 2k z2 d 2 = argr i2 + (2m + l)ir, 



(4) 



where m is an integer and arg z stands for the argument 
of the complex number z. Solving for d 2 , we get 
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argri2 - argr 23 + (2m + l)n 
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It is usually best to choose the value of m corresponding 
to the smallest positive admissible value of d 2 ; otherwise, 
internal resonances in the coating layer can spoil its an- 
tireflective properties for some angles of incidence. 

We shall now use condition (i) to determine the re- 
fractive index n 2 of the AR coating. Substituting the 
Fresnel's formulas (2) into the condition |r* 12 1 2 = |r* 2 3 1 2 , 
we obtain 

(S 2 - Sl ) 2 = (5 3 -5 2 )(5g-5 2 ) 
(Z 2 +Z 1 r (E 3 + E 2 )(E* 3 +E 2 y 
Straightforward algebra leads to 
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The s- and p-polarization cases need now to be consid- 
ered separately. Assuming materials 1 and 2 to be non- 
magnetic (fii = /i 2 = 1), in the s-polarization case we 
have 
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It is easily seen that Z 2 must be non-negative in order 
that k z2 be real, as we have assumed. For p polarization, 
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This leads to the quadratic equation for n 2 ,, 



n\ - Y%n\ + Y%n\ sin 2 (9 = 0, 



(11) 



which has real solutions 



ni = \ (V 2 2 ± ^/y 2 4 -4r 2 2 n?sin 2 ^ (12) 

provided that Y" 2 4 — 4Y 2 2 nf sin 2 # > 0. It can be shown 
that this condition, together with the condition of real- 
valuedness of k Z 2 [for both solutions of eq. (12)], is fulfilled 
if and only if 



Y" 2 2 > 4n? sin 2 9. 



(13) 



In practice, there are further constraints on the choice 
of the constituent material of the coating. Other exper- 
imental issues aside, rii is bounded from below by the 
refractive index of air, and from above, by the index at 
which a second propagative diffraction order appears at 
the given value of k x . In appendix A it is shown how 
such constraints of the general form 



2 ^ 2 2 
n min — n 1 — n max 



(14) 



can be transformed into equivalent constraints on the im- 
mittance of the PC, S3. Let us denote by S Q (a = 3, min, 
max) the reduced immittances X a /X\. It follows then 
that (14) holds if and only if one of the two following 
sets of conditions on S 3 is satisfied: 



£ 3 $ int P_ 
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or 
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S 3 £ int C m 



(15) 

where P- stands for the half-plane Re S 3 < 1 and C a 
(a = min, max) are the circles of radius 2 |1 — S 2 | cen- 
tered at ( \ ( 1 + S 2 ) , 0) . The symbols int A and ext A de- 
note the interior and exterior of a region A. To illustrate 
various possible geometries of the region of the complex 
S 3 plane determined by the constraints (14) transformed 
into the form (15), fig. 5 shows the shape of this region 
for s polarization and three distinct choices of the pa- 
rameters n m ; n , fimai, n>\ and 9. 

Step 2 The coating obtained in step 1 is not practi- 
cal, since its fabrication would call for integration of the 
PC with a completely different solid; moreover, a suit- 
able material with the required value of refractive index 
might not be easily available. However, as noted in sec- 
tion II, a homogeneous thin film can often be replaced 
without adverse effects by a subwavelength grating: this 
is step 2. Such a grating could be easily etched in the 
same process as the underlying PC; it would then natu- 
rally be composed of the same materials as the PC, with 
permittivities, say, c\ and £h (ei < £h)- 

In order to calculate the fill factor of a binary lamellar 
grating mimicking a layer with refractive index 712 ob- 
tained in the previous step, one can resort to the classical 



second-order effective-medium theory of gratings due to 
Rytov, described in ref. 14. According to this theory, the 
effective permittivity e s of a binary grating with period a 
and fill factor / (0 < / < 1), composed of materials with 
permittivities e\ and eh, and operating in the s polariza- 
tion is 
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where 



e s = fe h + (1 - /)ei. (17) 
For p polarization, the effective permittivity e p is 
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and e s is given by eq. (17). Thus, the required fill factor 
can be obtained by setting e s or e p to n 2 in eq. (16) or 
(18) and solving it numerically for /. 

It should be noted that in the domain of validity of 
Rytov's theory (small k^a/2'K) the functions e s (f) and 
e p (/) are monotonically increasing from e\ to eh- Thus, 
a binary grating cannot simulate a material with permit- 
tivity outside the range delimited by the permittivities 
of the grating's constituent materials. As a result, the 
bounds and n 2 ^ mentioned in step 1 must fulfill 
"■Lin > ei and n 2 max < e h , respectively. 

Step 3 The structure obtained at this stage should, 
in principle, ensure low reflectance for incidence an- 
gles close to 9. Nevertheless, owing to the applied 
approximations — neglect of higher diffraction orders ex- 
cited by the PC and the AR grating — its geometrical pa- 
rameters might not be precisely optimal. In addition, it is 
well known 13,14 that trapezoidal and triangular AR grat- 
ings have larger angular and frequency tolerance than 
lamellar ones. Therefore, as step 3, it is advisable to 
apply a numerical optimization procedure to adjust the 
geometry of the grating, described by some small number 
of parameters, so as to minimize a given objective func- 
tion p. The geometry obtained in step 2 can be expected 
to provide a good starting point for a local search algo- 
rithm, such as the Nelder-Mead simplex method (ref. 31, 
section 10.4). 



IV. APPLICATION TO A PHOTONIC 
CRYSTAL FLAT LENS 

The first PC we shall consider is a hexagonal lattice 
of air holes of radius r = 0.365a, where a is the lat- 
tice constant, etched in a dielectric matrix of permittivity 
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FIG. 5: Shaded: regions of the complex Z3 plane determined by the condition (15) equivalent to the constraint (14) for 
s polarization and (a) m = 1, n m i n = 1.5, n max = 3,9 — 0, (b) m = 1, n m in = 1, n. max = 3, 6 = 30°, (c) m = 1.5, n m i n = 1, 
n max = 3, 6 — 0. The circles C m i n and C ma x are defined in the text after eq. (15); note that in the case (b) C" max degenerates 
into the point (1,0). 




FIG. 6: Solid line: EFC of the PC considered in section IV at 
frequency u> = 0.311 x 2irc/a and for s polarization. Dashed 
line: EFC of air at the same frequency. 



e = 10.6. These parameters correspond to the structure 
whose fabrication was reported in ref. 32. For s polar- 
ization, at frequency uj = 0.311 x 2irc/a, the effective 
refractive index n = — 1 can be attributed to the crystal, 
since its equifrequency contour (EFC) takes an approxi- 
mately circular shape (fig. 6) with radius K « uj/c and 
group velocity directed inwards. (Information about the 
techniques used to obtain the numerical results discussed 
in this paper can be found in appendix B.) 

Veselago 37 predicted that a slab of material with n = 
— 1 should act as a flat lens: an image of an object placed 
near one of the surfaces of the slab should be produced on 
the other side of the slab. Figure 7 shows the map of the 
modulus of the electric field generated by a wire source 
with current 1 A located above a slab of the PC in ques- 
tion. The parts (a) and (b) refer to slabs truncated in 
the ways shown in figs. 8(a) and (b), respectively; from 
now on, these two structures will be referred to as SI 
and S2. In accordance with the theoretical predictions, 
images are formed below the slabs. However, their am- 
plitude is low (67 and 79V/m for structures SI and S2, 



respectively) and intense beams reflected from the top of 
the lenses are visible in the upper part of the plots. This 
suggests that only a small fraction of energy is transmit- 
ted through the lenses. Indeed, as shown in fig. 8(a), 
the reflectance of structure SI, \ro(0) | 2 , where ro is the 
specular reflection coefficient, exceeds 29% for all angles 
of incidence. Structure S2 performs better for low inci- 
dence angles, but degrades quickly with increasing 6. We 
shall now apply the algorithm presented in section III to 
design an AR grating for this PC, which ideally should 
function regardless of the angle of incidence. 

As discussed in refs. 38 and 39, the effective-medium 
model of PCs tends to be more accurate for crystals trun- 
cated along a plane with constant permittivity profile, 
as is the case, for instance, for structure SI. There- 
fore in the first step of the design procedure we shall 
calculate the effective transverse impedance Z^ of this 
structure. We consider two ways of obtaining this quan- 
tity. First, we calculate it in the framework of the model 
presented in ref. 38 using eq. (18) derived there. At 
frequency cu = 0.311 x 2iTc/a and angle of incidence 
9 = 45° (corresponding to k x = 0.220 x 2ir/a) we get 
Z x = 1.414 and Z 3 = 0.319. We should now check 
whether Z3 = Z%/Z\ = 0.225 lies within the region de- 
termined by the conditions (15) equivalent to the con- 
straints (14) with n min = 1 and n max = 2.51 (the max- 
imum index of a medium in which only a single prop- 
agative diffraction order exists). Figure 9, in which the 
value of Z-i cited above is marked with point A, shows 
that this is indeed the case. Therefore cqs. (5) and (9) 
can be used to calculate the parameters of the AR coating 
of the crystal: refractive index ni — 1.649 and thickness 
d 2 = 0.540a. The geometry of this structure, called S3 
from now on, is shown in fig. 10(a) and its reflectance 
is plotted in fig. 10(e) with a solid black line. It can be 
seen that the application of the coating reduces signifi- 
cantly the reflectance of the crystal, especially for small 
angles of incidence. However, the parameters of S3 are 
certainly not optimal, since its reflectance at the "design 
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FIG. 7: Modulus of the electric field generated by an s- 
polarized wire source with current 1 A located above a slab of 
the PC studied in section IV truncated along a plane (a) lying 
midway between two neighboring rows of holes, (b) crossing 
the centers of holes. 
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FIG. 8: Geometry of the PC studied in section IV truncated 
along a plane (a) lying midway between two neighboring rows 
of holes, (b) crossing the centers of holes, (c) Angular depen- 
dence of the reflectance of the structures shown in parts (a) 
and (b). 




FIG. 9: Shaded circle: region of the complex Z3 plane deter- 
mined by the condition (15) equivalent to the constraint (14) 
for s polarization, n\ = n m i n = 1 and n max = 2.51. Points A 
and B: reduced impedances Z3 of structure SI calculated in 
two different ways described in the text. 



angle" 9 = 45° is as large as 9%. This is due to the 
relatively large error introduced by the single-mode ap- 
proximation for negative-refraction PC bands, as pointed 
out in ref. 38. 

We shall evaluate, therefore, an alternative method of 
obtaining Z3, which consists in calculating it directly 
from the rigorous specular reflection coefficient r of 
the uncoated crystal at the chosen oj and k x . In other 
words, we assume that tq can be expressed in the form 
r = (Z 3 - Zi)/(Z 3 + Zi) [cf. eq. (2)] and invert this 
formula to obtain Z 3 = Zx(l + r )/(l — r ). Of course, 
the effective impedance defined in this way depends on 
the material properties of medium 1. Nevertheless, at 
least for the PC in question, this dependence is weak for 
sufficiently small m: we obtain Z 3 = 0.258 + 0.175i for 
n\ = 1 and the effective impedance does not change by 
more than 10% up to n\ = 2.25. As shown in fig. 9, the 
reduced impedance Z 3 — 0.182 + 0.124i corresponding to 
the above value of Z 3 (marked with point B) also lies 
within the allowed region of the Z 3 plane. Taking this 
value of Z 3l from eqs. (5) and (9) we get n 2 = 1.884 and 
d 2 = 0.565a. The angular dependence of the reflectance 
of the PC covered with this coating, shown in fig. 10(b) 
and called S4 in the following, is plotted in fig. 10(c) 
with a solid gray line. It is evident that this structure 
has much better angular tolerance than S3; moreover, 
its reflectance at 6 = 45° is only 0.05%. Therefore we 
choose S4 as a basis for the further steps of the algo- 
rithm. 

Numerical inversion of eq. (16) gives the fill factor / = 
0.192 of the binary grating mimicking a medium with 
n = 1.884. Since we would like the angular dependence 
of r to be symmetric with respect to 6 = 0, the grating 
should be positioned so as to preserve the vertical mirror 
symmetry axes of the underlying PC. This can be done in 
two ways, shown in figs. 10(c) and (d). The reflectance 
of these two structures, called S5 and S6, is plotted in 
fig. 10(e). Clearly, grating S5 reproduces fairly faithfully 
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FIG. 10: (a)-(b) Geometry of AR coatings S3 and S4, char- 
acterised by refractive index 712 and thickness specified 
next to the drawings, (c)-(d) Geometry of binary lamellar 
AR gratings S5 and S6, characterised by fill factor / and 
thickness e?2 specified next to the drawings, (e) Angular de- 
pendence of the reflectance of the structures shown in parts 
(a)-(d). 



the original reflectance curve of the AR coating S4. On 
the other hand, grating S6 behaves better in the high-# 
region. 

The lamellar gratings obtained in step 2 provide al- 
ready a remarkable improvement over the uncoated PC 
and, in contrast to the AR coatings from step 1, should 
be manufacturable. Nevertheless, their geometry can be 
further ameliorated. To this end, as mentioned in the last 
paragraph of section III, we use the Nelder-Mead simplex 
algorithm to find the optimum values of the dimensions 
u>i, w , hi, and h a parametrising the trapezoidal grat- 
ing shown in fig. 11. The objective function p is defined 
as the average of the numerically calculated reflectance 
of the given structure over the desired angular tolerance 
interval [^ m i„,S max ], 



P = 



^max ' 



f 

J 9rr 



d9. 



(20) 



The integral in eq. (20) is calculated with the 20-point 
Gauss-Legendre quadrature algorithm (ref. 31, section 
4.5), whose typical relative accuracy, ^10~ 5 , is better 




FIG. 11: Definition of the geometrical parameters Wi, w a , hi, 
and h of a trapezoidal grating superposed on the surface of 
structure SI. 



than that of the reflectance calculations, ^10~ 3 . The 
initial shape of the grating is taken to correspond to 
one of the lamellar gratings obtained in step 2, i.e., 
Wi = w = fa, hi = 0, and h Q = d 2 . The search rou- 
tine is terminated when the size of the simplex, defined 
as the average distance of its vertices from its geometric 
centre, falls below 10~ 5 . The final values of the geometri- 
cal parameters of the grating are determined by selecting 
the best among the 16 structures obtained by rounding 
each of the parameters delivered by the simplex algorithm 
upwards or downwards to a multiple of 0.01a. 

Application of this procedure with 9 m i n — 0° , 6> max = 
90° to grating S5 yields grating S7 shown in fig. 12(a). 
The plot in fig. 12(c) (solid black line) demonstrates the 
excellent antireflective properties of this structure (note 
the scale of the vertical axis). Its average reflectance is 
as low as 2.8%; in fact, |r o (0)| 2 does not exceed 5.5% 
until 9 = 87°. The structure does not seem to present 
special fabrication difficulties — e.g., acute angles — except 
possibly for the relatively thin dielectric veins separating 
the circular holes from the surface. Should this pose a 
real experimental difficulty, one can increase the value of 
hi at the expense of a slight performance deterioration. 
For example, grating S8 with hi — 0.08a [fig. 12(b)] has 
average reflectance of 4.8%. Figure 13 shows the map 
of the modulus of the electric field produced by a point 
source placed above a PC slab coated with AR gratings of 
type S7 from above and below. The comparison with fig. 
7 reveals the significant improvement brought about by 
the AR grating: not only are the reflected beams promi- 
nent in the upper part of the latter figure suppressed, but 
the amplitude of the image formed by the lens grows to 
159 V/m, which is two times better than in the situation 
from fig. 7(b). 

Optimization of structure S6 leads to gratings with av- 
erage reflectance comparable to that of S7 and S8 but 
composed of "narrower" trapezoids (wi+w ~ 0.2a), thus 
less suitable for fabrication. Therefore we omit the de- 
tailed discussion of these structures. 

As a complement, we mention that in refs. 40 and 
41 two other trapezoidal AR gratings, here denoted S9 
and S10, were presented. Their geometrical parameters, 
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FIG. 12: (a)-(d) Geometry of AR gratings S7-S10 character- 
ized by parameters Wi, w , hi, and h a specified next to the 
drawings, (e) Angular dependence of the reflectance of the 
structures shown in parts (a)-(d). To help visualize the de- 
tails of the |ro(#)| 2 dependence, the y axis has been truncated 
at I ml 2 = 0.2. 



shown in figs. 12(c) and (d), were obtained by minimizing 

the objective function | J^ 2 \r (6) \ dd (average modulus 
of the specular reflection coefficient r ) calculated with 
a less accurate quadrature algorithm. The average re- 
flectance of structure S9 in the full 0°-90° range, 2.8%, 
matches that of S7; in a more restricted range, say, 0°- 
80° , the performance of grating S9 is even slightly better. 
Nonetheless, its disadvantage lies in the presence of very 
thin dielectric veins at the surface. Grating SIO, with 
hi = 0.08a, is devoid of this problem. However, it is 
superseded by structure S8 with identical hi, which has 
somewhat lower average reflectance. 

Finally, a word about tolerance to fabrication imperfec- 
tions. A fabrication process invariably perturbs the ge- 
ometrical parameters of the manufactured structure. To 
assess the sensitivity of the proposed gratings to fabrica- 
tion errors, we have determined the maximum perturba- 
tion of each of the four geometrical parameters of grating 
S7 for which the grating's average reflectance in the angu- 
lar range < 9 < # max did not exceed 5%. Three values 
of #max were considered: 90° , 80° , and 60° . The results of 
this test are summarized in table I. It can be seen that the 
grating is more sensitive to variations of the height of the 
trapezoids (via the h a and hi parameters) than of their 
width (w and Wi). The constraints for # max = 90° and 
$max = 80° are rather stringent and unlikely to be met in 
practice. In contrast, fabrication of a structure satisfying 
the constraints for 9 max — 60° seems well within reach 
of current technology. We have also tested the frequency 
tolerance of grating S7, finding that the its average re- 
flectance stays below 5% for 0.3094 < uja/2nc < 0.3113 
(6> max = 90°), 0.3048 < u}a/2nc < 0.3122 (6» max = 80°), 
and 0.2511 < Loa/2-Kc < 0.3173 (6> max = 60°). This tol- 
erance seems quite sufficient for applications related to 
lensing. 
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FIG. 13: Modulus of the electric field generated by an s- 
polarized wire source with current 1 A located above a slab of 
the PC studied in section IV with S7-type gratings placed on 
its horizontal surfaces. 



V. APPLICATION TO A SUPERCOLLIMATING 
PHOTONIC CRYSTAL 

The second example to be considered is a PC composed 
of a square lattice of air holes of radius r — 0.3a, where 
a is the lattice constant, etched in a dielectric matrix of 
permittivity e = 12.25. Near the frequency oj — 0.265 x 
2-Kc/a its EFCs for p polarization take a square-like shape 
(cf. fig. 14). In consequence, supercollimated beams 42 ' 43 
can propagate in the crystal. 

The solid and dashed black curves in fig. 15(e) show 
the angular dependence of the reflectance of this PC at 
w = 0.265 x 27rc/ a and for two different truncation planes, 
corresponding to the structures shown in figs. 15(a)-(b) 
and called Sll and S12 in the following. It is seen that the 
crystal cut through hole centers has fairly low reflectance: 
about 10% at normal incidence and decreasing for larger 
angles up to 9 s=s 65°. This level of power losses might in 
fact be already sufficient for practical applications. Nev- 
ertheless, for the sake of illustration, we shall present the 
design procedure of AR gratings that help to decrease 
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flmax Wi Wo fa ho 

90 s 0.378-0.424a ( 22 nm) 0.151-0.170a ( 9nm) 0.045-0.055a ( 5nm) 0.543-0.557a ( 7nm) 
80° 0.330-0.463a ( 63nm) 0. 132-0. 186a (26 nm) 0.033-0.062a (14 nm) 0. 530-0. 569a (19 nm) 
60° 0. 192-0. 533a (162 nm) 0. 084-0. 208a (59 nm) -0.082a (39 nm) 0. 497-0. 586a (42 nm) 



TABLE I: Ranges of geometrical parameters of grating S7 for which its average reflectance at frequency 0.311 x 2-Kc/a in the 
angular range < 9 < 6> max does not exceed 5%. The numbers in parentheses are the lengths of the tolerance intervals for 
a = 476 nm, which corresponds to operation wavelength A = a/0.311 = 1530 nm. Note that the tolerance intervals correspond 
to perturbations of one parameter at a time (not all parameters simultaneously). 



0.50 




FIG. 14: EFC of the PC studied in section V at frequency 
uj — 0.265 x 2-Kc/a. The shaded region corresponds to the 
range 6 < 45° (\k x \ < 0.187 x 2n/a), where the EFC is ap- 
proximately flat and for which the minimization of the PC's 
reflectance is made. 



even further the reflectance of the PC in question. 

Here, we are mostly interested in coupling the incoming 
light to modes lying on the flat horizontal part of the PC's 
EFC. As shown in fig. 14, at frequency w = 0.265 x 2wc/a 
this corresponds roughly to the range \9\ < 45°, i.e., 
\k x \ < 0.187 x 27r/a. Therefore we choose 9 = 22.5° as 
the design angle of the AR coating. As in the previous 
section, we test two different ways of calculating the effec- 
tive immittance (in this case, admittance) of the crystal 
cut along a constant-permittivity plane, i.e., structure 
Sll. The effective-medium model presented in ref. 38 
yields Y3 = 6.138. In turn, the effective admittance cal- 
culated from the rigorous reflection coefficient of struc- 
ture Sll embedded in air is Y 3 = 6.075 — 1.191*. Figure 
16 shows that the reduced admittances corresponding to 
both these values lie within the region of the complex 
Y3 plane determined by the conditions (15) equivalent to 
the constraints (14) with n\ — n m ; n = 1, n max = 3.391 
(the refractive index for which the second propagative 
diffraction order appears) and 9 = 22.5°. The parame- 
ters of the AR coatings determined from these two val- 
ues of Y 3 are (n 2 = 2.548, d 2 = 0.374) and (n 2 = 2.595, 
d 2 = 0.391), respectively. Figures 15(c)-(d) show the ge- 
ometry of these coatings, henceforth referred to as S13 
and S14, whereas the angular dependence of their re- 
flectance is plotted with solid and dashed gray lines in 
fig. 15(e). As in the PC lens case, the AR coating S14 
designed using the value of Y3 obtained from the rigor- 
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FIG. 15: (a)-(b) Geometry of the PC studied in section V 
truncated along a plane (a) lying midway between two neigh- 
boring rows of holes, (b) crossing the centers of holes, (c)- 
(d) Geometry of AR coatings S13 and S14, characterized by 
refractive index 712 and thickness d2 specified next to the 
drawings, (e) Angular dependence of the reflectance of the 
structures shown in parts (a)-(d). 



ous reflection coefficient of the crystal performs slightly 
better than the other one. Therefore structure S14 shall 
be used in the subsequent design step. 

From numerical inversion of eq. (18) it follows that the 
fill factor of the binary grating mimicking a medium with 
n = 2.595 for p polarization is / = 0.812. Figures 17(a)- 
(b) show the geometry of the two gratings, called S15 and 
S16, with this fill factor and a mirror symmetry plane 
perpendicular to the direction of periodicity. From the 
juxtaposition of their reflectance curves [fig. 17(c)] it fol- 
lows that structure S15 has somewhat better performance 
than S16. In fact, there is some similarity between the 
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FIG. 16: Shaded circle: region of the complex Y 3 plane deter- 
mined by the condition (15) equivalent to the constraint (14) 
for p polarization, m = n m i n = 1 and n max = 3.391. Points 
A and B: reduced admittances Y3 of structure Sll calculated 
in two different ways described in the text. 



geometry of grating S15 and the truncated crystal S12, 
which also exhibited fairly low reflectance: the surface of 
both these structures contains "teeth" shifted by \a in 
the horizontal direction with respect to the positions of 
the circular holes. Therefore, one could view the crystal 
S12 as an imperfect realization of the AR grating S15. 

The lamellar grating S15 can be further ameliorated 
by adjusting its thickness d 2 and fill factor / to minimize 
the objective function p defined in cq. (20). We take 



= 0, 



as before, perform the opti- 



mization with the Nelder-Mead simplex algorithm. This 
leads to structure S17 with d 2 = 0.37a and / = 0.73, 
shown in fig. 17(c). Its reflectance curve is plotted in fig. 
17(d) (solid line). In the angular range < 9 < 45°, the 
reflectance never exceeds 0.6%, on average amounting to 
only 0.12%. The structure does not seem to present any 
special fabrication problems. It is possible to continue 
the grating's optimization by allowing it to take a trape- 
zoidal rather than a lamellar shape; however, in view of 
its already very good AR properties, this appears unnec- 
essary. 



VI. APPLICATION TO A NON-RECIPROCAL 
MIRROR 

In this section, the procedure is applied to a unidirec- 
tional mirror similar to that proposed in ref. 44. This 
device has the form of a slab of the PC shown in the in- 
set of fig. 18(a). The PC consists of a hexagonal lattice of 
non-centrosymmetric motifs etched in a magneto-optical 
matrix characterized by a gyrotropic permittivity tensor 
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FIG. 17: (a)-(c) Geometry of binary lamellar AR gratings 
S15, S16 and S17 characterized by fill factor / and thickness 62 
specified next to the drawings, (d) Angular dependence of the 
reflectance of the structures shown in parts (a)-(c) . 



the spatial inversion symmetry is also broken, the crystal 
becomes nonreciprocal: the usual property of the dis- 
persion relation, u>(k) = ui(—k) (ref. 45, pp. 22-23), no 
longer holds. The p-polarization EFC of this crystal at 
frequency lj — 0.4537 x 2-Kc/a is shown in fig. 18(a). 
It can be seen that plane waves impinging at the an- 
gle of about 44° (corresponding to k x — 0.315 x 27r/a) 
on the bottom surface of a slab made from this PC will 
be coupled to its propagative mode, and thus will be 
partially transmitted. However, waves travelling in the 
opposite direction, incident from the top, will be totally 
reflected, since the crystal has no propagative modes with 
-0.323 x 2?r/a < k x < -0.308 x 27r/a. This behaviour 
justifies the name unidirectional mirror. The problem 
with the presented device is its large forward loss: even 
waves propagating in the "allowed" direction undergo a 
significant reflection on the surfaces of the slab. This 
can be seen in fig. 18(b), where the reflectance curve of 
the crystal truncated as in fig. 18(a) is juxtaposed with 
the relevant fragment of the EFC. Clearly, to be useful 
in practice, the unidirectional mirror needs to be coated 
with some AR structure. 



The presence of this magneto-optical material lifts the 
time-reversal symmetry of Maxwell's equations. Since 



Unfortunately, the design of an appropriate AR grat- 
ing using the procedure described in section III turns out 
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FIG. 18: (a) p-polarization EFC of the magneto-optical PC 
shown in the inset at frequency cj = 0.4548 x 2nc/a. (b) Solid 
line: fc^-dependence of the reflectance |ro| 2 of this crystal, 
placed in air and truncated in the way indicated in the inset 
of part (a). Dashed line: a fragment of the EFC from part (a). 
The shaded region indicates the range of k x for which the 
crystal has a non-reciprocal gap, i.e. there are no propagative 
modes with the x component of the Bloch vector equal to — k x . 
(c) Shaded circle: region of the complex Yz plane determined 
by the condition (15) equivalent to the constraint (14) for 
p polarization, n\ — n m i n = 1 and n max = 1.50. The reduced 
admittance of the crystal, Ys = 13.5 — 7.3i, lies far beyond 
the range of the graph. 



to be impossible. For instance, at k x = 0.315 x 2irc/a, 
near the centre of the region where the unidirectional 
mirroring effect occurs, the reflection coefficient of the 
uncoated PC is tq = 0.900 — 0.055i, which corresponds 
to Y3 = 13.5 — 7.3z. As evidenced by fig. 18(d), this lies 
far outside the region of the complex Y3 plane determined 
by the conditions (15) equivalent to the constraints (14) 
with m = Timin = 1 arid ii ma x = 1-50 (the refractive index 
of the most optically dense coating in which the second 
propagative diffraction order would not yet appear). In 
fact, if we blindly apply eqs. (7) and (12) to calculate the 
refractive index of the optimum AR coating, we obtain 
rt2 = 5.79 or 0.70. It is obvious that none of these in- 
dices can be simulated by any binary grating composed 
of the constituent materials of the PC. Therefore, an AR 
structure for the unidirectional mirror will probably need 
to be designed with some purely numerical method. In 
particular, Lawrence et al. 2S have shown their approach 
to give good results for a superprism that, uncoated, has 
an extremely high reflectance (|?"o| 2 = 0.996). 



VII. EVANESCENT WAVES REGIME 

The behaviour of evanescent waves is approached 
through the example of the photonic crystal flat lens de- 
scribed in section IV. In that case, the amplification of 
these waves is of vital importance to obtain subwave- 
length resolution 2 . In a photonic crystal flat lens, two 
mechanisms may lead to an amplification of evanescent 
waves : the single interface resonances and the overall 
resonances (see ref. 46 for discussion in detail). While 
the first mechanism is actually involved in the original 
flat lens, the solely second mechanism has been investi- 
gated. This situation can be changed thanks to the new 
models introduced recently 28 ' 38,47 and based on reflectiv- 
ity on semi-infinite crystals. 

The present approach is based on the following con- 
jecture: assuming that a nearly constant effective index 
is obtained over all the range of wavevectors correspond- 
ing to propagative waves, it can be expected that this re- 
mains true when this range is continued in the evanescent 
regime. For instance, this behaviour has already been ob- 
served in multilayers 48 . For the photonic crystal flat lens, 
the refractive index is close to — 1 at ui = 0.311 x 2nc/a 
(from the dispersion law of Fig. 6) and, with the AR grat- 
ing, it can be considered that the effective index is nearly 
— 1 for the range of propagating waves. Now, if the effec- 
tive index remains nearly —1 for evanescent waves, then 
it can be expected that single-interface resonances are 
present around the working frequency oj = 0.311 x 2irc/a. 

Fig. 19 shows the single-interface resonances of struc- 
ture S7 below the light cone. It is observed that the 
dispersion curve of these resonances starts from the in- 
tersection between the bands of the bulk PC and the 
the light cone at cj = ck x = 0.311 x 2nc/a. This shows 
that the effective index has no jump from propagative 
to evanescent waves, and thus confirms the conjecture 
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FIG. 19: Single-interface resonances of the structure S7. The 
shaded areas are the bulk bands of the underlying PC. 



proposed above, as well as the result found in ref. 48 in 
the one-dimensional case. To our knowledge, no model 
can predict this continuity of the effective index at this 
change of regime. This might be a subject for furter in- 
vestigations. 



VIII. CONCLUSION 

In this article we have presented a new method of de- 
signing gratings that, superimposed on surfaces of PCs 
crystals, will minimize their reflectance. The design al- 
gorithm consists of three steps. First, the parameters of 
a homogeneous-layer AR coating are calculated from an 
effective-medium approximation of the PC in question. 
Second, an analytical effective-medium theory of gratings 
is used to find the parameters of a binary lamellar grat- 
ing composed solely of the constituent materials of the 
crystal and approximating the coating obtained in the 
previous step. Third, the shape of the grating is refined 
with a numerical local-search routine so as to minimize 
the reflectance of the structure in the desired angular or 
frequency range. This last step is necessary owing to the 
approximations made in the analytical derivations used 
in the first two steps of the procedure. 

This algorithm of AR grating design can be viewed as 
complementary to the method proposed by Lawrence et 
al. 27 ' 28 . While their approach is based on an exhaustive 
scan of the whole parameter space (made very efficient by 
the application of the matrix-valued effective immittancc 
of gratings), ours rests on approximate analytical consid- 
erations used to find a starting point for a local search 
procedure. 

The proposed method has been applied to three ex- 
ample crystals with EFCs of different curvature: a su- 
percollimating crystal with a very flat EFC, a crystal 
exhibiting negative refraction, with almost circular EFC, 
and a PC superprism, whose EFC has a kink. In the two 
first cases, the design process succeeded in producing AR 
gratings ensuring very low reflectance in a wide angular 
range. The obtained structures are quite compact and 



apparently rather straightforward to fabricate. In the 
last case, the procedure broke down owing to the viola- 
tion of the constraints (15) on the effective immittancc 
of the crystal that must be satisfied in order that the AR 
coating produced in the first step can be approximated 
with a binary grating made of realistic materials. The 
existence of these constraints is the basic limitation of 
the presented procedure. 



Appendix A: Constraints on H3 

We begin by noting that the conditions (14) can always 
be rewritten in the form 



~2 < ~2 < 772 

' — 'min — ' — 2 — "max 



(Al) 



with appropriate S m ; n and H max . Specifically, for s po- 
larization, the formulas for Z^ a and Z max follow readily 
from eq. (9): 



72 = 

min 



1 



C - n l sin2 1 



>7 £ _ 

max 2 



1 



(A2) 



n\ sin 



For p polarization, due to the more complex form of 
eq. (12) and the presence of the supplementary condi- 
tion (13), several cases must be considered. The final 
formulas for Y 2 - and Y 2 av are given in table II. 

nun nitix o 

To arrive at the form of the constraints on 23 , we sub- 
stitute eq. (7) into inequality (Al) and introduce reduced 
immittances S Q = S a /Hi (a = 3, min, max), obtaining 



S 2 . < 

mm — 



-ReS 3 + |~ 3 | 2 



ReS 3 -l 
This expression can be rewritten as 



(A3) 



Re Ha 



Re a, 



1 4- - 2 \ 

^P^J +(ImS 3 ) 2 
x(ReS 3 -l) > 0, 



1 + S 2 

1 1 — ' n 



(ImS 3 ) 2 - 



1 -H 2 



1 -S 2 



x(ReS 3 - 1) < 0. 

(A4) 

It follows that the constraints (14) are equivalent to the 
set of conditions written out in eq. (15). 



Appendix B: Numerical methods 

A word about the different methods used for the nu- 
merical simulations presented in this article. Almost 
all the calculations whose results are presented in sec- 
tions IV and VII have been made with the differen- 
tial method 33-36 . The exceptions are the field maps 



Range of n 2 



< nt < 



<C Tl x <C Tl m i n H max ^ 2t7- x 



n min n x ^ n 2 



mm 



"max 



< 9 mm x 4n^ 



.2 / „2 _2 < 2 2 



14 



TABLE II: Minimum and maximum bounds on Y 2 2 sufficient and necessary for fulfillment of the condition (14) together with 
the constraint (13) for at least one of the solutions (12) of eq. (11). The symbol n 2 denotes n\ sin 2 6. 



shown in figs. 7 and 13, produced with the finite- element 
method using the RF module of the COMSOL pack- 
age. The numerical results reported on in section V 
have been obtained with the finite-difference frequency- 
domain method. Finally, the simulations of the magneto- 



optical PC discussed in section VI have been made with 
a code based on the Fourier modal method 49 and imple- 
menting the factorisation rules for anisotropic materials 
derived in ref. 50. 
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